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Abstract 

In this paper, the numerical differentiation by integration method based on Jacobi polynomials originally introduced 
by Mboup, Fliess and Join (TjJ [2(| is revisited in the central case where the used integration window is centered. 
Such a method based on Jacobi polynomials was introduced through an algebraic approach (l9l |20| and extends the 
numerical differentiation by integration method introduced by Lanczos (1956) [2l[. The method proposed here, rooted 
in Til, 2(|, is used to estimate the nth (n G N) order derivative from noisy data of a smooth function belonging to at 
least C n+l+q (q G N). In [ljl,|2fj], where the causal and anti-causal cases were investigated, the mismodelling due to 



the truncation of the Taylor expansion was investigated and improved allowing a small time-delay in the derivative 
estimation. Here, for the central case, we show that the bias error is 0(h q+2 ) where h is the integration window length 
for / G (jn+q+2 - m |_ nc no j sc f rcc casc anc j fa e corresponding convergence rate is 0(S"+ 1 +i ) where 5 is the noise level for 
a well-chosen integration window length. Numerical examples show that this proposed method is stable and effective. 

Keywords: Numerical differentiation, Ill-pose problems, Jacobi orthogonal polynomials, Orthogonal series 



1. Introduction 



X 



Numerical differentiation is concerned with the numerical estimation of derivatives of an unknown function (defined 
from R to R) from its noisy measurement data. It has attracted a lot of attention from different points of view: observer 
design in the control literature EL 0, H i digital filter in signal processing 0, [B| , the Volterra integral equation of the first 
kind 0, 0] and identification . The problem of numerical differentiation is ill-posed in the sense that a small error 
in measurement data can induce a large error in the approximate derivatives. Therefore, various numerical methods 
have been developed to obtain stable algorithms more or less sensitive to additive noise. They mainly fall into five 
categories: the finite difference met ho ds 1 1 CI 1 1 ll . 1 1 2| . the mollification methods 13J, [lj, [15J , the regularization methods 



lil Il7l . |18| , the algebraic methods [19j, [20| that are the roots of the results reported here and the differentiation by 
integration methods [21, 22, 23 1, i.e. using the Lanczos generalized derivatives. 
The Lanczos generalized derivative Dhf, defined in [2l| by 



D h f(x) 



2h 3 



tf(x + t)dt 



-h 



3 

2h 



tf(x + ht) dt, 



is an approximation to the first derivative of / in the sense that Dhf(x 
of differentiation by integration. Rangarajana and al 



/' (x) + 0(h 2 ). It is aptly called a method 



22j generalized it for higher order derivatives with 



D h n) f{x) 



1 

h" 



j n L n (t)f(x + ht)dt, neN, 



where / is assumed to belong to C n+2 (I) with / being an open interval of R and L n is the nth order Legendrc 
polynomial. The coefficient 7„ is equal to lx3x5x ■^■x(2"+i) 2/1 > is the length of the integral window on which 
the estimates are calculated. By applying the scalar product of the Taylor expansion of / at x with L n they showed 
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that D^ ] f{x) = f^(x) + 0{h 2 ). Recently, by using Richardson extrapolation Wang and al. [23| have improved the 
convergence rate for obtaining high order Lanczos derivatives with the following afhne schemes for any neN 



(t) (o„ f(x + ht) + b n f(x + X n ht)) dt, 




where / is assumed to belong to C" +4 (I), a n , b n and A„ are chosen such that f(x) = f^ n \x) + 0(h 4 

Very recently an algebraic setting for numerical differentiation of noisy signals was introduced in [24I ] and analyzed 



The reader may find additional theoretical foundations in [251 . 1 261 ] . The al geb raic manipulations used in 
are inspired by the one used in the algebraic parametric estimation techniques [27 , H, 2{j . Let us recall that 
|20H analyze a causal and an anti-causal version of numerical differentiation by integration method based on Jacobi 
polynomials 

where / is assumed to belong to C n (I) with / being an open interval of R. The coefficient 7 K)M) „ is equal to 
(-1)" 

(m+")K™+")! ' wnere are tw0 integer parameters and h > is the length of the integral window on which 
the estimates are calculated. In [l9| the authors show that the mismodelling due to the truncation of the Taylor 
expansion is improved allowing small time-delay in the derivative estimation. Here in this article, we propose to 
extend these differentiation by integration methods by using as in 1|| 2^ Jacobi polynomials: for this we use a central 



estimator (the integration window is now [—1, 1]) and the design parameters are now allowed to be reals which are 
strictly greater than —1. It is worth to mention that in most of the practical applications the noise can be seen as an 
integrable bounded function (which noise level is 8 as is considered in this paper). Another point of view concerning 
the noise definition/characterization is given in j2f| for which unbounded noise may appear. Let us mention that 
the Legendre polynomials arc one particular class of Jacobi polynomials, that were used in [22[ and [23| to obtain 
higher order derivative estimations. Moreover, it can be seen that these so obtained derivative estimators correspond 
to truncated terms in the Jacobi orthogonal series. In fact, the choice of the Jacobi polynomials comes from algebraic 
manipulations introduced in the recent papers by M. Mboup, C. Join and M. Fliess [ljl [2(|, where the derivative 
estimations were given by some parameters in the causal and anti-causal cases. Here, we give the derivative estima- 
tions in the central case with the same but extended parameters used in Til 2(3|. If / G C' n+q+2 then we show that 



the bias error is 0(h q+2 ) in the noise free case (where 2h is the integration window length). We also show that the 

q+l 

corresponding convergence rate is 0(<5™+ 1 +«) for a well-chosen integration window length in the noisy case, where S 
is the noise level. One can see that the obtained causal estimators in [f9|, H3] are well suited for on-line estimation 
(which is of importance in signal processing, automatic control, etc.) whereas here the proposed central estimators 
are only suited for off-line applications. Let us emphasize that those techniques exhibit good robustness properties 
with respect to corrupting noises (see [2^, [3(| for more theoretical details) . These robustness properties have already 
been confirmed by numerous computer simulations and several laboratory experiments. Hence, the robustness of the 
derivative estimators presented in this paper can be ensured as shown by the results and simulations reported here. 

This paper is organized as follows: in Section [2] firstly a family of central estimators of the derivatives for higher 
orders are introduced by using the nth order Jacobi polynomials. The corresponding conver gen ce rate is 0(h) and 
can be improved to 0(h 2 ) when the Jacobi polynomials are ultraspherical polynomials (see [31(). Secondly, a new 
family of estimators arc given. They can be written as an affine combination of the estimators proposed previously. 
Consequently, we show that if / G C n+1+q (I) with geN the corresponding convergence rate is improved to 0(h q+1 ). 
Moreover, when the Jacobi polynomials are ultraspherical polynomials, if / G C n+2+q (I) for any even integer q the 
corresponding convergence rate can be improved to 0(h q+2 ). Numerical tests are given in Section [3] to verify the 
efficiency and the stability of the proposed estimators. 



2. Derivative estimations by using Jacobi orthogonal series 

Let f 5 = f + zu be a noisy function defined in an open interval I C R, where / G C n+1 (I) with n G N and the 

noiscQ zu is bounded and integrable with a noise level 5, i.e. 5 = sup |tu(:e)|. Contrary to [22| where the nth order 

xei 



1 More generally, the noise is a stochastic process, which is bounded with certain probability and integrable in the sense of convergence 
in mean square. 



Legendre polynomials were used, we propose to use, as in [19|, |20|, the nth order Jacobi polynomial so as to obtain 
estimates of the nth order derivative of /. The nth order Jacobi polynomials (see 3l|) are defined as follows 



3=0 



n + a\ ( n + 0\ ft — 1 
n-j 



t + 1 



(1) 



where a, (3 e]-l,+oo[. Let us denote Vc/i , g 2 € C°([-l,l]), (gi(-) 7 g 2 (-)) a /3 = j\w a ^(t)g 1 (t)g 2 (t)dt, where w aJ3 (t) = 
(1 — t) a (l + t) 13 is the weight function. Hence, we can denote its associated norm by || • \\ a ,/3- 

We assume in this article that the parameter h > and we denote Ih '■= {x £ I, [x — h,x + h] C /}. 

2.1. Minimal estimators 

In this subsection, let us ignore the noise va for a moment. Then we can define a family of central estimators of 

/( „) 

Proposition 2.1 Let f 6 C n+1 (I), then a family of central estimators of f( n ' can be given as follows 

V* G h, D%l p f(x) = ^J Pn, a At) f( x + M)dt, (2) 

where p n ^(t) = ^ZlZZti) ^^ ™ th B(n + a + l,n + (3 + I) = E&S+g^gBl . 

Moreover, we have D h n ^ a af(x) = f^ n \x) + 0(h). 

Remark 1 In order to compute p n , a .p, we should calculate P^"' whose computational complexity is 0(n 2 ). Hence, 
the computational effort of p n ,a,i3 "is 0(n 2 ). 



Proof. By taking the Taylor expansion of /, we obtain for any x £ Ih that there exists 9 €]x — h, x + h[ such that 



Un-in Ln+ljn+l 

f(x + ht) = f{x) + htfix) + ■■■ + -^f {n \x) + /<" +1 >(0). 

n! (n + lj! 



(3) 



Substituting §3§ in ([5]), we deduce from the classical orthogonal properties of the Jacobi polynomials (see |3l[) that 

r.1 



Pn,a,p(t) t m dt = 0, < m < n, 



Pn^At) t n dt = in\). 



(4) 
(5) 



Using ([3]) and ([5]), we can conclude that 

= / <wW /(* + ht)dt = + o(h). 

Hence, this proof is completed. 



□ 



In fact, we have taken an nth order truncation in the Taylor expansion of / in Proposition 12 . 1 1 where n is the order 
of the estimated derivative. Thus, we call these estimators minimal estimators (see 3, 2fJ]). Then, we can deduce 
the following corollary. 

Corollary 2.2 Let f S C n+1 (I), then by assuming that there exists M n +i > such that for any x £ I, (x) | < 

Af n +i, we have 



D%l, f(x)-r>(x) 



(«)/ 



<Cyh, 



(6) 



wfcere Ci = |i n+1 p«,a^(i)| dt. 



When a = j3 the Jacobi polynomials are called ultraspherical polynomials (see 31|). In this case, we can improve 



Lemma 2.3 Let Pn"'^ be the nth order ultraspherical polynomials, then we have 



1 Pi a ' a \t)w a , a (t)t n+l dt = 0, (7) 
-l 



where I is an odd integer. 

Proof. Recall the Rodrigues formula (see (3l| ) 



P^\t)w a At) = ^r^K+^+n(*)i. ( 8 ) 

we get, by substituting ((5J in ([7]) and applying n times integrations by parts, that 

i (*)«;„,,?(«) <ft = £ «, a+rii/3+n (t) i< dfc. (9) 

If a = /3 and Z is an odd number then w a+ri: p +n (t) t l is an odd function and the integral in ^ is equal to zero. Hence, 
this proof is completed. □ 



Consequently, we can deduce from Proposition 12.11 and Lemma 12.31 the following corollary. 
Corollary 2.4 Let f € C n+2 (I) and a = (3 in Provosition \2.1[ then we obtain 

e h, d%Ij(x) = + o(h*). (io) 

Moreover, if we assume that there exists M n+ 2 > such that for any x 6 /, |/ < -™ +2 '( x )| ^ -Wn+2j we have 

D { :lj{x)-f {n \x) <C x h\ (11) 

OO 

where C 1 = |i"+ 2 pn,a,a(0l ^' 

We can see in the following proposition the relation between minimal estimators of / and minimal estimators of 
/ (n) . 

Proposition 2.5 Let f G C n+1 {I), then we have 

*^ ^-sl^S^C) (l2) 

where a n j = a + n — j and /3,- = j3 + j. 

In order to prove this proposition, we give the following lemma. 

Lemma 2.6 for any i £ N, we /lave 



(P^\t),f{x + ht)) i /\ 9 -, 

VXG7 ' 1 ' || P (^)||2 "2^ ^ U , . (1 . ; . i U M - 



where Qfjj = a + i — j and /3j = /3 + j. 

Proof. Observe from the expression of the Jacobi polynomials given in ([T]) that 



we get 



(p^\t),f(x + ht)) a E °) (*. + fj J\ ^.ft (*) /(* + &*) (15) 

Then, by using Proposition 12.11 with n — and Pq Qi,j ' ,3j ' > (t) e 1 we obtain 



Recall that (see |3l[) 

Hp(a,/3),|2 2°+*+ 1 r(q + z + l)r(/3 + z + l) 

11 1 2i + a + /3 + ir(a + /3 + i + l)r(i + l)' 1 ' 

then the proof is completed by using (jTT)) in (|T5|) . □ 

Proof of Proposition [2751 From it is easy to show after some calculations that 

n W 1 r(q + /3 + 2n+l) \ v ^ v V a ^ 

"h,a,pTK x ) (2h) n T{a + l3 + n+l) \\Pk a,f>) f a ' 



Hence, this proof can be completed by using Lemma 12761 and (fTojl . □ 



Now, we can see in the following proposition that the estimates given in Proposition 12.11 are also equal to the first 
term in the Jacobi orthogonal series expansion of fW(x + ht) at point t = 0. 

Proposition 2.7 Let f £ C n+1 (I), then the minimal estimators of given in Proposition \2. 1\ can be also written 
as follows 

(p( a+n 'P +n \t),fW(x + ht)\ 
V* e I h , = ± ||pW+n)| , '-^^ Pt +n '" +n) (0). (19) 

Moreover, we have 

V* e I h , D^jfix) = D^ +n ^ +n f^(x). (20) 
Proof. By using the Rodrigues formula in ^ and applying n times integrations by parts we get 

^ (_])n,2-(2n+a+0+l) pi Jra 
Drl R f(x) = j r / — \w a+n B+n(t)} f(x + ht) dt 



h n B(n + a + 1, n + /3 + 1) J_ x df 

2-(2n+Q+^+l) /■! 



B(n + a + l,n + /3 + l) y_j 
n(°) f( n )M 



U7 a +n^+n(t)/ (n) (a; + At)di 



Then, by using p( a + n >P+ n ) (t) = 1 and \\P^ +n ' 0+n) \\l +n ^ +n = 2 2n+a +^ +1 B{n + a + 1, n + + 1), we can achieve 
this proof. 

□ 

2.2. Affine estimators 

It is shown in Proposition ^ . 71 that the minimal estimators of f( n \x) given in Proposition ^ . 1 1 are equal to the value 
of the order truncated Jacobi orthogonal series expansion of f( n \x + ht) at t = 0. Let us assume that / € C ,l+1 (I), 
then we define now the gth (q £ N) order truncated Jacobi orthogonal series of f( n \x + ht) by the following operator 

9 (pt^'^^J^Hx + h-)) 

^ e A, + := E ^"^M- (21) 



Take t = in (|2ip. we obtain a family of estimators of f( n '(x) with 



(a+ra,/9+ra) 112 



+n,f)+n p(a+n,/3+n) 



(0). 



(22) 



i=0 II- 1 i lla+n,/3+n 

To better explain our method, let us recall some well-known facts. We consider the subspace of C°([— 1, 1]), defined 



by 



U — n,> Q „ / p(a+n./3+n) p (a+n^+n) p (a+n,/3+n) \ 

rtq — span < _r , r x , ■ ■ ■ , r q > 



Equipped with the inner product (•, -) a+n p +n , H q is clearly a reproducing kernel Hilbert space [32 
reproducing kernel 

1 p(a+n,/J+n) / \-p(a+n,fl+7i) t 



I p (a+n,P+n) |, 2 

The reproducing property implies that for any function f^ n '(x + h-) belonging to C°([— 1, 1]), we have 



(lC q (;t),fM(x + h-)) =D%lpJ(x + th), 

\ I a+7iJ3-\-n 



(23) 
with the 

(24) 



(25) 



where D^ a - f(x + h-) stands for the orthogonal projection of f^ n \x + h-) on H q . Thus, the estimators given in ([22 



can be obtained by taking t = 0. 

We will see in the following proposition that the estimators Djf^ a q f(x) can be written as an afHne combination 
of different minimal estimators. These estimators are called affine estimators as in [lij . 

Proposition 2.8 Let f G C n+1 (I), then we have 



(o+n,/3+n) , .2i + a + (3 + 2n + 1 



'(0) 



t + a + 8 + 2n + 



where q £ N, a^j = a + i — j cmd /3j = (3 + j . Moreover, we have 

~* D (a+n,/3+n) , 2i + a + ft + 2n + I 



'(0) 



i + a + f3 + 2n + 



3=0 



™ =1- 



(27) 



Proof. By replacing a by a + n, (3 by /3 + n and /(a; + ht) by f^ n '{x + ht) in Lemma T2.61 we obtain 



LP. 



(a+n,/3+n) ,, 2 



3=0 



i\ 2i + a + (3 + 2n + 1 



f (28) 



j/ i + a + /3 + 2n + 



II Q+n,a+n 

Then (|26|) can be obtained by using (|20|) and (|22j) . By using the Binomial relation, (|27[) can be easily obtained. □ 



Hence, by using Proposition 12.11 an explicit formulation of these afhne estimators is obtained in the following 
corollary. 



Corollary 2.9 Let / € C n+1 (L), then the affine estimators of can be written as 



V* G h, D%lpJ(x) = I Q a ,p, n ,S) /(* + ht ) dt, 



where 



(29) 



(30) 



with n S N, 0„ , a, ozwen in Proposition \2.1\ and on = a + i — 7, 8j = + i. 



Consequently these afnne estimators are also differentiation by integration estimators. 



Remark 2 Q a ,f),n,q * s a sum of ^(q + l)(q + 2) terms. According to Remark\^ the computational effort of each term 
is 0{n 2 ). Hence, the computational effort of Q a ,/3,n,q is also 0(n 2 ). 

It is shown in Proposition 12 . II that the convergence rate of minimal estimators is 0(h). We will see in the following 
proposition that the convergence rate of affinc estimators can be improved to 0(h q+1 ). 



Proposition 2.10 Let f G C n+1+q (I) with q G N, then we have 

G h, D%lpJ(x) = /<»)(*) + 0(h q+1 ). 



(31) 



Moreover, if we assume that there exists M n+ i +q > such that for any x G /, |/'™ +1+IJ, )(x)| < AI„+i+ q , then we have 



in), 



(32) 



where C * = ^T^fl I-i \t n+1+q Q a ^, hq (t)\dt. 

Proof. By taking the Taylor expansion of /, we get for any x G Ih there exists £ g]x — h, x + h[ such that 

un+l+q+n+l+q 



f(x + ht) = f n+q (x + ht) + 



(n + ! + <?)! 



(n+l+q) 



(0, 



(33) 



n+q 



-\ h J P 

where f n + q (x + ht) = — ~f (x) is the (n + q)th order truncated Taylor series expansion of fix + ht). 



3=0 



Let us take the Jacobi orthogonal series expansion of f^ q (x + ht). Then by taking t — 0, we obtain 

f(«) frA -Y^J I / a+n,l3+n 



i=0 



j(a+n,/3+rt) i|2 



I a-\-n,a.-\-n 



Similarly to (|22|) we obtain 



fn+q( x ) ~~ ~frn j Qa,P,n,q(t)fn+q( 



ht)dt, 



(34) 



(35) 



from (|2^|) where Q a ,p,n, q is given in Corollary 12.91 by ([50)1 . 

By calculating the value of the nth order derivative of f^ q at t = 0, we obtain f%+ q (x) = f^ n \x). Then by using 
(1291) and (ESI) we obtain 



D h!L,pJ(*) - f (n) ( X ) =w I . [/(* + ht) - f n+q (x + ht)} dt 



1 

h q+1 



(n + l + q)\ J_ t 
-~0{h q+1 ). 



Consequently, if for any x £ I |/( n + 1 +9) (x) | < Mrj+i+g, then we have 

r-l 



D ( :i,j(x)-f (n \x) 



< h q - 



(n + l + g)! 7_ 



t n+1+9 g a) ^„, g (t)l^. 



□ 



We can deduce that the affine estimator for f^ n \x) obtained by taking the gth order truncated Jacobi orthogonal 
series expansion of f( n '(x + h-) can be also obtained by taking the (n + q)th order truncated Taylor series expansion 
of / with a scalar product of Jacobi polynomials. 



n+q 



hPV 



Moreover, let f n+q (x+ht) = f n (x+ht)+r q (x+ht) where r q (x+ht) = — — f^\ x ) for 1 _t 1 an d r q (x+ht) = 

j=n+l J ' 

for g = 0, then becomes 



9 1 

where P = 

i=o 



/&(*)=£ 

2=0 

(a+n,/3+n) / ,\ (n) 



IP, 



(a+n,^+n) ,i 2 



+n,p+n p(a+n,/3+n) /qn _|_ ^ 



(t),ri" ; (x + W) 



IP 



(a+n^+n) |i 2 



(0). 



Observe that fl^ } {x -\-h-) is a 0th order polynomial, then by using the orthogonal properties of P i 



(a+n,/3+n) 



we have 



i (P^ P+n) (t),ft\x + ht)) 



i=0 



IP 



(a+n,P+n) 112 



(0) 



a-\-n,a-\-n 

(t),ti n >(x + ht) 



I p(a+n,£+n) i|2 
I II Q+n,a+n 



By calculating the value of the nth order derivative of f^ q and /™"^ at t = 0, we obtain f£+ q (x) = fn (x) = (x). 
Hence, we get R = 0. Hence, we can deduce that 



P 



1 

17' 



Qa,p,n,q{t)r q (x + ht)dt = 0, 



(36) 



where Q a ,p, n ,q is given in Corollary 12.91 by ([50]) . 

Consequently, (|36p explains why the convergence rate can be improved from 0(h) to 0(/i 9+1 ): the price to pay is 
some more smoothness hypotheses on the function /. 

If we consider the noisy function f s , then it is sufficient to replace f(x + h-) in ([29)1 by f s (x + h-) so as to estimate 
fW(x). Then we have the following definition. 

Definition 2.11 Let f s = f + vo be a noisy function, where f € C n+1 (J) and w is a bounded and integrable noise 
with a noise level 5. Then a family of estimators of f^ is defined as follows 



Vxel h> Di n l p q f 5 (x) = ^ I Q a , , n , q (t)f s (x + ht)dt, 



(37) 



where Q a ,p,n,q is given by (jffflp . 

In the following proposition we study the estimation error for these estimators. 

Proposition 2.12 Let f s be a noisy function where f G C n+1+q (I) and vo is a bounded and integrable noise with a 
noise level 5, then 



D 



i:Uqf d (x)-f w (x) 



(»)/ 



< C 2 h" +1 + c 3 



(38) 



where Ci is given in Proposition ^. l0\ and C3 = |_j \Q a .p,n.q(t)\ dt. 
Moreover, if we choose h 



"C 3 x 
(q+l)C 2 ° 



then we have 



p(n) 



Proof. Since 



\dF> a J S (x) - D[ n l a J(x)\\ = \\d ( A a _ \f s (x) - f(x)] II < — f 1 |Q„ „ „ Jt)\ dt 



(39) 



by using Proposition 12 . 101 we get 



D 



x < D^J\x)-D^ q f(x) D^J(x)-f^(x) 

where C3 = |Q a „3,n,g(*)| Let us denote the error bound by ip(h) = C2h q+1 + C^-^. Consequently, we can 

1 

calculate its minimum value. It is obtained for h* = , rrfc 5 5 and 



fS 



(n) 



,(«) 



_(q+l)C 2 

n +l + q fq+1' 
9+1 



C 2 " +1+q C 3 " +1+? 



Then, the proof is completed. 



(40) 

□ 



In Proposition 12. 8[ we improve the convergence rate from O(h) to 0(h q+1 ) (q G N) for the exact function / by 
taking an affinc combination of minimal estimators of f( n \ Here, the convergence rate is also improved for noisy 

1 q+l 1 

functions. It passes from 0(S"+ I ) to 0(S n + 1 +i) if we choose h = cS n + 1 +i , where c is a constant. 

Remark 3 Usually, the sampling data are given in discrete case. We should use a numerical integration method to 
approximate the integrals in our estimators. This numerical method will produce a numerical error. Hence, we always 
set the value of h larger than the optimal one calculated in the previous proof. 

We have seen in the previous subsection that the convergence rate of minimal estimators can be improved to 0{h 2 ) 
when a = (3. Let us then study the convergence rate of affinc estimators in this case. 



Corollary 2.13 Let f G C n+2+q (I) where q is an even integer. If we set a = f3 in i2fy) , then we have 



(41) 



Moreover, if we assume that there exists M n+ 2+ q > such that for any x G /. | /^ n+9+2 ^ (a;) | < M n+2 +q, then we have 



(«), 



< C 2 h q+2 , 



(42) 



where & = (t+g+2)! I-i \t n+q+2 Q^n, q {t)\dt and 



Q««w = E4 a+n ' ft+n) (o)E(-D 



2i\ 4z + 2a + 2n + 1 
j ) 2i + 2a + 2n + 1 



LP; 



I a+n,a+n 



(43) 



i=0 j=0 

with p n .a 2i j,fij given in Proposition \2.1\ and a 2 i,j = a + 2i — j , f3j = a + j . 

Proof. Observe that P^ n ' a+n) (-t) = (-l)(«+ x ) P £+ n ' a + n ) ( t ) f or any 4 e [_i f 1] ( see [HI p.80), we obtain 

p(a+n, a+n )^ = q HcncCi ^ becomes 

r)( n ) f It) — / a+w,a+ra p (a+n,a+n) ,„x 

^h^a.qJ W n n (a+n,Q+n)|| 2 



If / 6 C n+2+9 (I), then let us take f n + q +i as the (n + q + l)th order truncated Taylor series expansion of f(x + ht). 



r (n) 



By taking the Jacobi orthogonal series expansion of f n+q+ \ 



q+l ( P (^, a+n){t)J (n^ +i{x + ht) 



i=0 



(a+n,a+n) 1 1 2 



(0), 



we obtain 



p(n) 



(n) 



P (n) 



ll a+n,a+n) i|2 ' 



IP: 



I a+n,a+n 



/ Qa 1 a,n 1 g(<)[/(a + /lt)-/n+ g +l(x + fa)]dt 



-1 



(n + 2 + g )! 
0(/i 9+2 ). 



y Q a , a , n , q (t)t n+2+q f {n+2+q) {?)<&, C e]x -h,x + h[ 



Consequently, (|42|) follows directly from the hypothesis on |/(™ +9+2 )(x)| . Since pj a+n ' a+n > (jfj = o for any odd integer 
i, (|43l) can be obtained by using (|30p. Then this proof is completed. □ 

Remark 4 According to fs2 l. we can deduce the asymptotic behavior of the number £' when h — > + 

\e-x\ _ 



lim 



1 



h^o+ h n + q + 3 
Similarly to Proposition 12 . 1 2l we can obtain the following corollary. 



(44) 



Corollary 2.14 Let f G C n+2+q (I) where q is an even integer. If a = /3 in Definition \2.11[ then the estimation error 
f or D h? a , a ,qf 5 ( x ) is 9 iven b y 

where Ci is given in Corollary \2.13\ and C3 is given in Proposition \2.12\ 
Moreover, if we choose h 



nC, 



(q+2)C 2 



th 



en we nave 



00 

In the following proposition, if we assume that / G C ,l_1 (I) then we can define the generalized derivative of f( n \ 
We can see that if the right and left hand derivatives for the nth order exist, then this generalized derivative converges 
to the average value of these one-sided derivatives. 

Proposition 2.15 Let f 6 C n_1 (J), then we define the generalized derivative of by 



V.t 6 h, D%i a J(x) = / Q a , n , q (t) f(x + ht) dt 



where Q a ,n.q is defined by Moreover, if /j. (x) and /^(x) exist at any point x G Ih, then we have 



(n), 



(45) 



(46) 



where /l" (resp. f_) denotes the right (resp. left) hand derivative for the nth order. 
Before proving this proposition, let us give the following lemma. 

Lemma 2.16 Let n G N and Q a ,n.q be the function defined on [—1, 1] by \43§ where q is an even integer. If n is even 
then Q a .n,q is also even, odd else. 



Proof. By taking a = j3 in (|22p . we obtain 



p 



(a+n,a+n) 



(0) 



— ' II p( Q +"- Q +")||2 
=0 W^i \\a+n,a+n 



Pt +n ' a+n) (t)w a+n , a+n (t)fM(x + ht)dt. 



By using (|14[) and replacing a, /3 by a + n, we get for Z = 0, • • • , n — 1 



(47) 



P 



(a+n,a+n) 



(-2)' 



E 



i + a + n \ / z + a + n 



t-3 



where ai+nj = a + i + n. — j, o*j+ n = a + j + n. Then, by applying the Rodrigues formula, we get 



dt l 

if 

3=0 



t 


- a - 




(i + a + n\ 




J 




V i-3 ) 



where a>i+ n -i.j = a + i + n — j — I, aj+ n -i = a + j + n — l. Hence, we get that 4p [P i 
at —1 and 1. Thus, by applying n times integrations by parts in (|47[) . we obtain 



d' rp(a+n,a+n) 



w a+ n. a +n] are equal to 



P 



(a+n,a+n) 



(o) r 1 



By using Corollary 12.91 with a = /5, we get 



(a+n,a+n) n 2 

II ot-\-n ,oi-\-n 



Pi («*+»,aH-n) (t)tu _^ (t) f( X + ht)dt. (48) 



qw*) = (-i)"E 



p(a+n,a+n) ^ ^ n 
(a+n,a+n)||2 



Pt +n ' a+n) (t)w a+n , a+n (t) 



Since P i ( - a+n ' Q+ '^(o) = for any odd integer i, (|4"9"|) becomes 



(49) 



Qa,n,,(*) = (-l) n X] 



-j(a+n,a+ro) 
2i 



(0) d" 



P 



(a-\-n,a-\-n) 



2i 



{t)w a + n ,a+n(t) 



Since p^ +n ' a+n ^_^ = p(a+n,a+n) ^ ^- gee j^jj-j an( j l y a+n)a+n (_t) = io a+n)Ce+n (t), we have 



<7f' 



P 2 \ ; (i)w a+ „, a+n (t) =(-1) — P 2 \ ; (-i)Wa+n,a+ri(-0 



Thus, we have Q a ,n,q{t) = (— l) n <9a,n,g( — *)• Then this proof is completed. 



(50) 



□ 



Proof of Proposition I2TT51 Let us recall the local Taylor formula with the Peano remainder term [35(. For any 

given e' > 0, there exists S > such that 



f(x + ht)-f n -i(x + ht) 



f^(x) 



(hty 



and 



fix + ht) - f n _ x { x + ht) - i+±l(hty 



< e'\ht\ n , for S < ht < 0, 



< e'(ht) n , for < ht < 6, 



(51) 



(52) 



where f n _i(x + ht) is the (n — l)th order truncated Taylor series expansion of f(x + ht). Let us consider the function 
g(x) = x n the nth order derivative of which is equal to (n\). Thus, by using (|22[) we have 

V* G 4, D^n(x) = (n!). 



Thus, by applying Corollary 12.91 with a = (3, we get 



1 f 

V ^ e 4, £»lS, ai ^(a:) = Q a ,n >q {t)9{x + = (n!). 

In particular, by taking x = we get f\Q a ,n, q (t) (ht) n dt = (71!). According to Lemma [2.161 t n Q ajl ^ q {t) with 
£ G [—1, 1] is an odd function. Hence, we have J_ 1 Q a ,n,q(t) (ht) n dt = ^ J Q a ,n,q(t) (ht) n dt. Thus, we get 



ga, n , 9 (t)^-^(^)"dt = ^/i' l) (*), 
! n! 2 



and 



By using ([22)) and Corollary 12.91 with a — (H we get 

1 f 1 

Va- G I h , D^l a q f n ^(x) = jjij Q a ,n, q (t) U-i{x + ht) dt = 0. 
Hence, by using (|53|) . (|54]l and (|55|) we obtain 

^ fl /(*)-|(/i n) (*)+/^(*)) 



<- 



Qa,n,q(t) f(x + ht) - f n -l(x + ht) 



f H (x) 



(htr 



,, q (t) \ f(x + ht) - f n _x(x + ht) - i±-±-L(ht)' 



dt 



dt. 



By using (|43"|) . we get 



\Q*,n„(t)\dt<±\i£^*\o)\£(*y 



3=0 



j j 2i + 2a + 2n + 1 



Pn,a 2i::i ,0j (t) \ dt. 



(53) 



(54) 



(55) 



(56) 



(57) 



Then, according to © and (fl"4")l we can obtain that L \p n , a2 i (*) | ^ < 00 ■ Hence, 

/ |Qa,n, fl (*)t n |d*< / |Qa,n, 9 (*)| <&<«>. 

Jo Jo 

Consequently, for any e > 0, by using (|56|) . (|5T|) and ([52l with £ = 2e' J |Q Q ,n,g(i) i n | d£j there exists <5 such that 
< h < 6 and 



< £. 



Then, this proof can be completed. 



□ 



3. Numerical tests 



In order to demonstrate the efficiency and the stability of the previously proposed estimators, we present some 
numerical results in this section. First of all, we analyze the choice of parameters for these estimators. 




(a) n = 1, q = 0,2, • • • , 8 and a = 0, 1, • • • , 10. (b) n = 2, g = 0, 2, • • • , 8 and a = 0, 1, • • • ,10. 



Figure 1: Values of Cz- 



'°a,o C 4 lo a,o c 3 

3r , , 4.5 r 




(a) q = 4, n = 1, • • • , 4 and a = 0, 1, • • • , 10. (b) q = 4, n = 1, • • • ,4 and a = 0, 1, • • • ,10. 



Figure 2: Values of log 10 C4 and log 10 C3. 



3.1. Analysis for parameters' choice for the bias term error and the noise error 

As is shown previously, the proposed estimators contain two sources of errors: the bias term error which is produced 
by the truncation of the Jacobi orthogonal series expansion and the noise error contribution. The error bounds for 
these errors are given in Corollary |2.14l We are going to use the knowledge of the parameters' influence to these error 
bounds. This will help us to obtain a tendency on the influence of these parameters on the estimation errors. 

According to Corollary 12.131 we set a = f3 and choose the truncation order q to be an even integer. On the 
one hand, it is clear that we should set q as large as possible so as to improve the convergence rate and reduce the 
bias term error. On the other hand, the noise error contribution is bounded in Corollary 12.141 by B no i se = C^-j^ 

where C3 = J_ t \Q a . n ,q{t) \ dt. We can see in Figure [T] the different values of C3 where n = 1, 2, q = 0, 2, • ■ • , 8 and 
a = 0, 1, • • • , 10. It is clear that with the same values for n and a, C3 increases with respect to q. Furthermore, it is 
easy to verify that C3 increases with respect to q, independently of n and a. Hence, in order to reduce the bias term 
error and to avoid a large noise error, we set q = 4 in our estimators. With this value, according to Corollary 12 . 141 the 
convergence rate is 0{8 ™+ 6 ). 

The bias term error is bounded by B bias = C 2 h^+ 2 in Corollary E33 where C* 2 = ™+g+ 2 y. j\ \t n+q+2 Q a ^ q {t)\ dt. 

Let us introduce C 4 = J\ \t n+q+2 Q a , n , q (t)\dt. We can see in Figure [2] the different values of log 10 C4 and log 10 C3 
when n = 1, • • ■ ,4, q = 4 and a = 0, 1, ■ • • ,10. It is clear that C4 decreases with respect to a while C3 increases with 
resnect to a. Thus, in order to reduce the bia,s term error, we should set n as larse a,s nossible. However, a, larse value 



of a may produce a large noise error contribution. Here, we choose a = 5. 

Until here, we have chosen q = 4 and a = 5. The noise error decreases with respect to h and the bias term error 
increases with respect to h. In the next subsection we are going to choose an appropriate value for h by using the 
knowledge of function / and by taking into account the numerical integration method error. 




(c) n = 3, a = 5, q = 4 and h = 777T 3 . (d) n = 4, a = 5, q = 4 and h = 850T S . 



Figure 3: The exact values of f^ n (xt) and the estimated values Dj^' a a q fi(xi) for <5 = 0.15. 
3.2. Simulation results 

The tests are performed by using Matlab R2007b. Let f s (xi) = f(xi) + cw{xi) be a generated noise data with an 
equidistant sampling period T s = where c > 0. The noise cw{x.i) arc simulated from a zero-mean white Gaussian 
iid sequence by the Matlab function 'randn' with STATE reset to 0. By using the well-known three-sigma rule, we 
can assume that the noise level for cm is equal to 3c. We use the trapezoidal method to approximate the integrals in 
our estimators with 2m + 1 values. The estimated derivatives of / at the point x% & I = [—2, 2] are calculated from 
the noise data f s (xj) with Xj £ [— x, — h,X{ + h], where h = mT s and 2m + 1 is the number of sampling data used 
to calculate our estimation inside the sliding integration windows. When all the parameters are chosen, Q a ,/3, n ,q in 
the integrals of our estimators can be calculated explicitly by off-line work with the 0(n 2 ) complexity. Hence, our 
estimators can be written like a discrete convolution product of these pre-calculated coefficients. Thus, we only need 
2m + 1 multiplications and 2m additions to calculate each estimation. 

The numerical integration method has an approximation error. Thus, the total error for our estimators can be 
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(a) n = 1, a = 5, q = 4 and h = 442T S 



(b) n = 2, a = 5, q = 4 and h = 549T S 
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(c) n = 3, a = 5, q = 4 and h = 643T S . 



(d) n = 4, a = 5, q = 4 and h = 733T S 



Figure 4: The exact values of f^™' (i;) and the estimated values a q f2(%i) for <5 = 0.15. 



bounded by 

T m (QWO f\xi + k-)) - / (n) (Xi)| < |T m (Q«, n , g (.) / 5 fa, + h-)) - T m (Q a , n , q (-) /fa + h-))\ 

+ \T n (Q„,n,,(-) /fa + ft')) - £j!la, g /(3i)| + \ D kLJ&) - /^fa) 
^ Bnoise ~t~ B num -\- Bijias ^total 3 

where T m (Q a . n ^ q (-) f{x t + h-)) (rcsp. T m (Q ai „ )g (-) /'(a* + ft-))) is the numerical approximation to D^ ^f 'fa) 

(resp. q g /lfa)) with the trapezoidal method and B num is the well-known error bound for the numerical 

integration error |36| : 

^Hc^/fa) - T ™ (Qa,«,,(-) /fa + M) 



< 



10 / x 2 Sup (Q at n,q(t) /fa + ^)) (2) = B„ 

12(2m)^ t e [_i,i] 



(58) 



We are going to set the value of m such that B to tai reaches its minimum and consequently the total errors in the 
following two examples can be minimized. For this, we need to calculate some values of fW*> with k = 0, • ■ ■ , n + q + 2. 



According to Rcmark|4j we calculate the value of M n +2+q in the interval 
the function / is unknown. 



n+g+3 5 n+q+3i 



However, in practice, 



Example 1. We choose fi(x) = sin(27ra:)e x as the exact function. The numerical results are shown in Figure [3j 
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(a) n = 1, a = 5, q = 4 and h = 442T S 
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(b) n = 2, a = 5, q = 4 and h = 549T S 
Errors for ^ 
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(c) n = 3, a = 5, q = 4 and h = 643T S . 



(d) n = 4, a = 5, <j = 4 and h = 733T S 



Figure 5: The estimation errors for the estimated values a /2(^i) for 8 = 0.15 
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(a) The estimation errors for D £ n ^ 5 4 /2(a;i) with varying val- 
ues of hi for <5 = 0.15. 



(b) Values of m; 



Figure 6: Errors for improved estimations with varying values of hi. 
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(a) Noise error with its error bounds 



(b) Bias term error + numerical integration error with the error 
bounds 



Figure 7: The estimation errors and their corresponding error bounds for D ' s 5 4/2(2)1) with varying values of hi for 8 = 0.15 



and the dash-dotted lines represent the estimated derivative values D^f^ a g fi(xi). Moreover, we give in Table Q] the 
total error values max ^ D^ a a q f\{xi) — f\ { x i) f° r the following noise levels: 5 = 0.15 and 5 = 0.015. We can see 
also the total error values produced with a larger sampling period T' s = 10T S = 10~ 2 . 



Table 1: max D« ( Xi ) - /<»> ( Xi ) 



5 


n = 1 (to) 


n = 2 (m) 


n = 3 (m) 


n = 4 (m) 


0.15 


9.45e-002 (591) 


1.1 (698) 


1.258e + 001 (777) 


1.278e + 002 (850) 


0.015 


1.85e - 002 (425) 


2.951e-001 (523) 


3.888 (601) 


4.588e + 001 (675) 


0.015 (Tj = 0.01) 


4.06e - 002 (47) 


5.645e - 001 (55) 


7.359 (62) 


9.686e + 001 (69) 



Example 2. When fi{x) = e x , we give our numerical results in FigureU]with the noise level S = 0.15, where the cor- 



for 5 = 0.15 and 6 = 0.015, where the total error values are produced with T s and a larger sampling period T' = 10 -2 



responding errors are given in Figure [5] In Table [21 we also give the total error values max D^' a fii^i) — f% ( x i) 

Xi£[2,2] '■ ' ' q 



Table 2: max D " J 2 { Xi ) - f^ n> (x t ) 

Xi G [2,2] I ' ' ' 



s 


n = 1 (m) 


n = 2 (to) 


n = 3 (to) 


n = 4 (to) 


0.15 


1.42e-001 (442) 


2.152 (549) 


2.982e + 001 (643) 


3.756e + 002 (733) 


0.015 


2.22e - 002 (346) 


4.435e - 001 (428) 


5.973 (510) 


8.769e + 001 (595) 


0.015 (Tj = 0.01) 


3.404e - 001 (54) 


3.425 (61) 


3.638e + 001 (68) 


5.235e + 002 (79) 



We can see in Figure [5] that the maximum of the total error for each estimation (solid line) is produced nearby 
the extremities where the bias term error plus the numerical error (dash line) are much larger than the noise error. 
The noise error (dash-dotted line) is much larger elsewhere. This is due to the fact that the total error bound B to tai 
is calculated globally in the interval \—h — 2, 2 + h]. The value of m with which B to tai reaches its minimum is used 
for all the estimations D^ a a f2{xi) with Xi £ [—2, 2]. This value is only appropriate for the estimations nearby the 
extremities, but not for the others. In fact, when the bias term error and the numerical integration error decrease, we 
should increase the value of m so as to reduce the noise errors. In order to improve our estimations, we can choose 
locally the value of m = nii, i.e. we search the value rrn which minimizes B to tai on [—hi + Xi, Xi + hi] where hi = rriiT s . 
We can see in Figure |H] the errors for these improved estimations £)£ ™ 5 5 4 f{xi). The different values of TOj are also 
given in Figure [6] The corresponding error bounds are given in Figure [7] We can observe that the error bounds 
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(a) 5 = 0.15, n = 1, a = 5, q = 4 and h = 1700T S 



(b) S = 0.015, n = 1, a = 5, q = 4 and h = 1200T S 





(c) 5 = 0.15, n = 2, a = 5, q = 4 and h = 1700T S . 



(d) <5 = 0.015, n = 2, o = 5, q = 4 and h = 1200T S 
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(c) 8 = 0.15, n = 3, a = 2, q = 2 and h = 1700T S 



(f) <5 = 0.015, n = 3, a = 2, g = 2 and h = 1500T S 



Figure 8: The exact values of f^"\xi) and the estimated values D Q q f3(xi)- 



help us to know the tendency of errors so as to choose parameters for our estimations. On the one hand, the chosen 
parameters may not be optimal, but as we have seen in our examples, they give good estimations. On the other hand, 
the optimal parameters q op , a op and m op with which the total error bound reaches its minimum may not give the best 

i^oti T~n a f i r\n T'ticit io ^rtnr wr* nnlu men ttincn nrvnr Kmmrlc tn plinnco tlio TT-ciln^ nm 



Example 3. Let us consider the following function 



--U 3 + 2a;, ifx<0, 

/a(*) = < f , 

-af 5 + 2x, if .t > 0, 
6 



which is C 2 on / = [—2, 2]. The second derivative of jfe is equal to \x\. Consequently, docs not exist at x = 0. If 
n > 1, then this function does not satisfy the condition / e C n+2+q (I) of Corollary 12.141 The numerical results are 
shown in FigurelH where the sampling period is T s = 10~ 3 and the noise level 5 is equal to 0.15 and 0.015 respectively. 
The solid lines represent the exact derivative values of /|"^ for n = 1,2,3 and the dash-dotted lines represent the 
estimated derivative values D^ a a q fz{xi). For the estimations of and f( 2 \ we set a — 5 and q = 4. When wc 
estimate f^ 3 \ the noise error increases. Hence, we need to decrease the values of a and q to a = 2 and q = 2. In Tabic 



[31 we give also the total error values max 

XiG[2,2] 



00 

h,a,a,q 



for 



1,2 and 6 = 0.015,0.15. 



Table 3: max l-D^^/aCxi) - / 3 (n) (x») 



5 


n = 1 (m) 


n — 2 (to) 


0.15 


9.7e- 003 (1700) 


9.65e - 002 (1700) 


0.015 


4.7e - 003 (1200) 


7.23e- 002 (1200) 
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